###--------------------------------------------------------------------###
#
# R Codes for the practicals of "Fundamentals in Ecology"
# Visualization on R 2022
# Solution to the tasks
# 
# Bachelor 2 SIE - Sciences et Ingenierie de l'Environnement
# Ecole polytechnique federale de Lausanne (EPFL)
#
# by Eugenie Mas, Christoph Bachofen, Hannes Peter, Jonas Paccolat
#
###--------------------------------------------------------------------###

### Tasks 1 Debugging:
### Change the code so that it works. What has to be changed? Type in the correct code that works. 

# # 1. using basic plot, plot the length of the sepal from the iris dataset in function of the petal length from the the dataset iris
#plot(iris$Petal.Length, iris$Sepal.Length) 

# # 2. using ggplot2, do the same plot
ggplot(iris, aes(x=Petal.Length, y=Sepal.Length))+
geom_point()

# # 3. now you want to color by Species
ggplot(iris, aes(x=Petal.Length, y=Sepal.Length,color=Species))+
geom_point()

# # 4. actually you would like to see only the sepal length between 0 and 6cm
ggplot(iris, aes(x=Petal.Length, y=Sepal.Length))+
geom_point()+
ylim(0,6)

# # 5. you want to change also the name of the y axis 
ggplot(iris, aes(x=Petal.Length, y=Sepal.Length))+
geom_point()+
scale_y_continuous(name= "Sepal length (cm)", limits=c(0,6))


###--------------------------------------------------------------------###

### Task 2:
#You are a scientist working on the effect of cold temperature on the performance of the grass species Echinochloa crus-galli.
#You collected data of CO2 uptake by the plants in function of the ambient CO2 of plants growing under cold temperatures and ambient temperatures and providing from Qu?bec or Mississippi.
#At the end of your experiment, you end up in a dataframe of 84 rows and 5 columns that you can find pre-loaded in R under the name CO2.
#It contains the following objects:

#Plant: an ordered factor with levels Qn1 < Qn2 < Qn3 < ... < Mc1 giving a unique identifier for each plant.
#Type: a factor with levels Quebec Mississippi giving the origin of the plant
#Treatment: a factor with levels nonchilled chilled
#conc: a numeric vector of ambient carbon dioxide concentrations (mL/L).
#uptake: a numeric vector of carbon dioxide uptake rates (umol/m^2 sec).

#1) First you want to see if the CO2 uptake if impacted by the treatment chilled/nonchilled, plot the most accurate plot.
#Q) what conclusion can you do? 
#answer: The plants assimilate generally less under colder temperature. 
library(ggplot2)
ggplot(CO2, aes(x=Treatment,y=uptake))+
  geom_boxplot()

#2) However, your plants comes from different area, you want to see if the origin of the plant influences the effect observed. Plot the most accurate plot.

#Q: what conclusion can you do? 
#Answer: the Qu?bec plants seem more adapted to cold temperature, as the decrease in CO2 assimilation seems smaller for the Qu?bec plants.
ggplot(CO2, aes(x=Treatment,y=uptake,color=Type))+
  geom_boxplot()

#Q: Is there some outliers that seems to appear?
#Answer: Yes in Quebec plants, the lowest values seem to be outliers.

### Task 3: 
#You decide to present the concentration uptake with a barplot.

#1) Do a barplot of the CO2 uptake in function of the treatment and the origin of the plant
ggplot(CO2, aes(x=Treatment, y=uptake,fill=Type))+
  geom_bar(stat="identity", position=position_dodge())

#1) This is great but actually it would be nicer to add some error bar. 

#First you have to calculate the sd and mean per treatment and type
#It exist several options to Calculate the average and standard deviation of the uptake CO2 per condition of temperature and the origin of the plants
#and to obtain in the end a dataframe with 4 lines with in each line there is one treatment, one origin and the corresponding mean and sd
#for example, the first line should be "non chilled / Quebec / 35 / 9"
#Here I present 4 options

#Option 1: The most DIY one
Quebec<-subset(CO2,Type=="Quebec") #split the dataset per Type
Mississippi<-subset(CO2,Type="Mississippi")

Quebecmean<-aggregate(uptake ~ Treatment, data=Quebec, mean)#calculate the mean of treatment for Quebec
Quebecsd<-c(by(Quebec$uptake, Quebec$Treatment,sd)) #calculate the sd of treatment for Quebec
barplotQuebec<-cbind(Quebecmean,Quebecsd) #combine the vector in a dataframe
colnames(barplotQuebec)<-c("Treatment","mean","sd") #rename the column just to be more clear
barplotQuebec$Type<-c("Quebec","Quebec")#add a column with the type

#do the same manip but with the mississippi dataset
Mississippimean<-aggregate(uptake ~ Treatment, data=Mississippi, mean)#calculate the mean of treatment for Mississippi
Mississippisd<-c(by(Mississippi$uptake, Mississippi$Treatment,sd)) #calculate the sd of treatment for Mississippi
barplotMississippi<-cbind(Mississippimean,Mississippisd) #combine the vector in a dataframe
colnames(barplotMississippi)<-c("Treatment","mean","sd") #rename the column just to be more clear
barplotMississippi$Type<-c("Mississippi","Mississippi")#add a column with the type

#combine the 2 datasets mississippibarplot and quebec barplot
barplot<-rbind(barplotQuebec,barplotMississippi) #combine the line of 2 datasets 


#Option 2: Less DIY but still long
CO2$code<-paste(CO2$Treatment, CO2$Type) #create a unique code that is the treatment + the Type in a new column (to have in this column 4 different factors)
co2mean<-aggregate(uptake ~ code, data=CO2, mean)#calculate the mean using this new code
co2sd<-c(by(CO2$uptake, CO2$code,sd)) #calculate the sd using this new code
barplot2<-cbind(co2mean,co2sd) #combine the vector in a dataframe
colnames(barplot2)<-c("code","mean","sd") #rename the column just to be more clear

# to reproduce the same barplot than in question 1, you have to add a column treatment and a column type
#there is 2 options to do that: option a: 
barplot2$Treatment<-c("chilled","chilled","nonchilled","nonchilled") 
barplot2$Type<-c("Mississippi","Quebec","Mississippi","Quebec")

#option b: the prettiest: you just split the "code" column with the function separate from the package tidyr
library(tidyr)
barplot2<-barplot2 %>% separate(code, c('Treatment', 'Type'))


#Option 3: the more efficient using "aggregate"

co2mean2<-aggregate(uptake ~ Treatment+Type, data=CO2, mean) #calculate the mean of treatment per type
co2sd2<-aggregate(uptake ~ Treatment+Type, data=CO2, sd)#calculate the sd of treatment per type
barplot3<-cbind(co2mean2,co2sd2) #combine the two dataframe
barplot3<-barplot3[,-c(4,5)] #remove the column the repeat the info, just to make it prettier
colnames(barplot3)<-c("Treatment","Type","mean","sd") #rename the columns to not mess up the columns

#option 4: the lazy solution but the most efficient one: use the function "summarySE" in the packge Rmisc, it automatically calculate the mean and sd for the treatment and the origin of the plant
library(Rmisc)
barplot4<-summarySE(CO2,measurevar="uptake",groupvars=c("Treatment","Type"),na.rm=T)

#Finally, you can create your plot with ggplot2:

ggplot(barplot, aes(x=Treatment, y=mean, fill=Type))+
  geom_bar(stat="identity",position=position_dodge())+
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2,
                position=position_dodge(.9),size=1)


### Task 4:
#Actually, you are more interested in the relationship between the the CO2 uptake and the ambiant CO2 concentration

#1) PLot this relationship
#Q: What is the tendency?
#Answer: The Co2 uptake by the plant increase with the ambient co2 concentration.
ggplot(CO2, aes(x=uptake,y=conc))+
  geom_point()

#2) Highlight in this plot the treatments with different colors
ggplot(CO2, aes(x=uptake,y=conc,color=Treatment))+
  geom_point()

#3) Add a linear regression. Do you see a difference between treatments?
#Answer: the relationship is less strong for the chilled treatment
ggplot(CO2, aes(x=uptake,y=conc,color=Treatment))+
  geom_point()+
  geom_smooth(method="lm", fullrange=T,se=F, size=2)

#4) You want to see if the origin of the plant impacts this relationship, split your plot in function of the origin. Is it change something?
#answer: the Quebec plants seems to react less strongly to the treatment than the Mississippi plants
ggplot(CO2, aes(x=uptake,y=conc,color=Treatment))+
  geom_point()+
  geom_smooth(method="lm", fullrange=T,se=F, size=2)+
  facet_wrap(~Type)

#Bonus question: Let's go to the next level and upgrade your plot to do the most beautiful plot that you can and post it on the moodle, let see who did the prettiest! :)

ggplot(CO2, aes(x=uptake, y=conc, color=Treatment))+
  geom_point(size=2)+
  scale_color_manual(values=c("darkorange","darkcyan"))+
  scale_x_continuous(name="Ambiant CO2 concentration (mL/L)") +
  scale_y_continuous(name="CO2 uptake rates (umol/m? sec)")+
  geom_smooth(method="lm", fullrange=T,se=F)+
  theme_test(base_size = 30)+ #the option "base_size" allows to increase everything in the plot: the symbols, name, letter etc
  facet_wrap(~ Type)